---
title: Statistical analysis
subtitle: Litter decomposition
author: <a href="marceloarayasalas.weebly.com/">Marcelo Araya-Salas</a> & Andrea Vincent
date: "`r Sys.Date()`"
toc: true
toc-depth: 2
toc-location: left
number-sections: true
highlight-style: pygments
format:
html:
df-print: kable
code-fold: true
code-tools: true
css: qmd.css
editor_options:
chunk_output_type: console
---
```{=html}
<style>
body
{ counter-reset: source-line 0; }
pre.numberSource code
{ counter-reset: none; }
</style>
```
```{r set root directory, echo = FALSE}
# set working directory as project directory or one directory above,
rootdir <- try(rprojroot::find_rstudio_root_file(), silent = TRUE)
if (is(rootdir, "try-error")) rootdir <- ".."
knitr::opts_knit$set(root.dir = rootdir)
```
```{r add link to github repo, echo = FALSE, results='asis'}
# print link to github repo if any
if (file.exists("./.git/config")){
config <- readLines("./.git/config")
url <- grep("url", config, value = TRUE)
url <- gsub("\\turl = |.git$", "", url)
cat("\nSource code and data found at [", url, "](", url, ")", sep = "")
}
```
```{r setup style, echo = FALSE, message = FALSE, warning=FALSE}
# options to customize chunk outputs
knitr::opts_chunk$set(
class.source = "numberLines lineAnchors", # for code line numbers
tidy.opts = list(width.cutoff = 65),
tidy = TRUE,
message = FALSE,
warning = FALSE
)
```
<!-- skyblue box -->
<div class="alert alert-info">
# Purpose {.unnumbered .unlisted}
- Evaluate role of nutrient availability on litter decomposition
</div>
```{r load packages, echo = FALSE, message = FALSE, warning=FALSE}
# remotes::install_github("rlesur/klippy")
# github packages must include user name ("user/package")
pkgs <- c("kableExtra", "knitr", "readxl", "ggplot2", "tidybayes", "cowplot", "ggrepel", "posterior", "ggridges", "viridis", "brms", "ggdist")
# install/ load packages
sketchy::load_packages(pkgs, quite = TRUE, upgrade.deps = TRUE)
```
```{r functions and global parameters, echo = TRUE, message = FALSE, warning=FALSE}
cols <- viridis(10, alpha = 0.7)
fill_color <- viridis(10)[7]
# brms models
chains <- 4
iters <- 40000
prior <- c(prior(normal(0, 10), "b"), prior(normal(0, 50), "Intercept"),
prior(student_t(3, 0, 20), "sd"))
# set ggplot2 them
ggplot2::theme_set(theme_classic(base_size = 20))
# standard error
se <- function(x) sd(x) / sqrt(length(x))
```
```{r, echo=FALSE}
source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
source("~/Dropbox/R_package_testing/brmsish/R/check_rds_fits.R")
```
# Nutrient change
## Nitrogen
Content
```{r}
nutr <- read.csv("./data/raw/litter-nutrients-mas.csv")
nutr$plot.f <- as.factor(nutr$plot)
nutr$days.sc <- scale(nutr$days)
nutr$litter.n.content.prop.initial <- nutr$litter.n.content.perc.initial / 100
# remove plot 9
sub.nutr <- nutr[nutr$plot != 9, ]
agg_n <- aggregate(litter.n.content.perc.initial ~ colecta + treat + days, nutr, mean)
agg_n$sd <- aggregate(litter.n.content.perc.initial ~ colecta + treat + days, nutr, sd)$litter.n.content.perc.initial
agg_n$se <- aggregate(litter.n.content.perc.initial ~ colecta + treat + days, nutr, se)$litter.n.content.perc.initial
agg_n$treat <- factor(agg_n$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_n, aes(x = days, y = litter.n.content.perc.initial, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = litter.n.content.perc.initial + se, ymin = litter.n.content.perc.initial - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter N content (% initial)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_n$days),
labels = unique(agg_n$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/Litter_N_content.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/Litter_N_content.tiff)
```{r, eval = FALSE}
fit.n <- brm(litter.n.content.prop.initial ~ treat * days.sc + (1 | plot.f), data = nutr, chains = chains, cores = chains, family = Beta(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), prior = prior, file = "./data/processed/litter.n.content.rem_model", file_refit = "on_change")
```
```{r, results='asis', warning=FALSE}
extended_summary(read.file = "./data/processed/litter.n.content.rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
Concentration
```{r}
agg_conc_n <- aggregate(perc.n ~ colecta + treat + days, nutr, mean)
agg_conc_n$sd <- aggregate(perc.n ~ colecta + treat + days, nutr, sd)$perc.n
agg_conc_n$se <- aggregate(perc.n ~ colecta + treat + days, nutr, se)$perc.n
agg_conc_n$treat <- factor(agg_conc_n$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_conc_n, aes(x = days, y = perc.n, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc.n + se, ymin = perc.n - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter N concentration (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_conc_n$days),
labels = unique(agg_conc_n$days)) + theme(legend.position = c(0.3, 0.8))
ggsave("./output/Litter_N_concentration.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/Litter_N_concentration.tiff)
```{r, eval = FALSE}
nutr$prop.n <- nutr$perc.n / 100
fit.perc.n <- brm(prop.n ~ treat * days.sc + (1 | plot.f), data = nutr, chains = chains, core = chains, family = Beta(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), prior = prior, file = "./data/processed/litter.n.perc_model", file_refit = "on_change")
```
```{r, results='asis'}
extended_summary(read.file = "./data/processed/litter.n.perc_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
## Phosphorus
Content
```{r}
nutr$litter.p.content.prop.initial <- nutr$litter.p.content.perc.initial / 100
# excluding plot 9
agg_p <- aggregate(litter.p.content.perc.initial ~ colecta + treat + days, sub.nutr, mean)
agg_p$sd <- aggregate(litter.p.content.perc.initial ~ colecta + treat + days, sub.nutr, sd)$litter.p.content.perc.initial
agg_p$se <- aggregate(litter.p.content.perc.initial ~ colecta + treat + days, sub.nutr, se)$litter.p.content.perc.initial
agg_p$treat <- factor(agg_p$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_p, aes(x = days, y = litter.p.content.perc.initial, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = litter.p.content.perc.initial + se, ymin = litter.p.content.perc.initial - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter P content (% initial)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_p$days),
labels = unique(agg_p$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/Litter_P_content.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/Litter_P_content.tiff)
```{r, eval = FALSE}
fit.p <- brm(litter.p.content.prop.initial ~ treat * days.sc + (1 | plot.f), data = nutr, chains = chains, prior = prior, cores = chains, iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), file = "./data/processed/litter.p.content.rem_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/litter.p.content.rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
Concentration
```{r}
agg_conc_p <- aggregate(perc.p ~ colecta + treat + days, sub.nutr, mean)
agg_conc_p$sd <- aggregate(perc.p ~ colecta + treat + days, sub.nutr, sd)$perc.p
agg_conc_p$se <- aggregate(perc.p ~ colecta + treat + days, sub.nutr, se)$perc.p
agg_conc_p$treat <- factor(agg_conc_p$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_conc_p, aes(x = days, y = perc.p, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc.p + se, ymin = perc.p - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter P concentration (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_conc_p$days),
labels = unique(agg_conc_p$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/Litter_P_concentration.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/Litter_P_concentration.tiff)
```{r, eval = FALSE}
# convert to proportions to use beta distribution
sub.nutr$prop.p <- sub.nutr$perc.p / 100
fit.perc.p <- brm(prop.p ~ treat * days.sc + (1 | plot.f), data = sub.nutr, chains = chains, cores = chains, family = Beta(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), prior = prior, file = "./data/processed/litter.p.perc_model", file_refit = "on_change")
```
```{r, results='asis'}
extended_summary(read.file = "./data/processed/litter.p.perc_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
<div class="alert alert-success">
# Takeaways {.unnumbered .unlisted}
- Litter P content is significantly lower in plus N treatment plots than in control plot, after accounting for variation explained by time
</div>
## Remaining litter
```{r, eval = TRUE, out.width = "100%", echo = FALSE, fig.align= "center"}
dat <- read_excel("./data/raw/litter_data.xlsx")
dat$plot.f <- as.factor(dat$plot)
dat$prop.litter.rem <- dat$perc.litter.rem / 100
dat$days.sc <- scale(dat$days)
agg_rem <- aggregate(perc.litter.rem ~ colecta + trat + days, dat, mean)
agg_rem$sd <- aggregate(perc.litter.rem ~ colecta + trat + days, dat, sd)$perc.litter.rem
agg_rem$se <- aggregate(perc.litter.rem ~ colecta + trat + days, dat, se)$perc.litter.rem
agg_rem <- agg_rem[order(agg_rem$trat), ]
pd <- position_dodge(15)
ggplot(agg_rem, aes(x = days, y = perc.litter.rem, color = trat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc.litter.rem + se, ymin = perc.litter.rem - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter mass remaining (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_rem$days),
labels = unique(agg_rem$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/Litter_mass_remaining.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/Litter_mass_remaining.tiff)
```{r, eval = FALSE}
fit <- brm(prop.litter.rem ~ trat * days.sc + (1 | plot.f), data = dat, chains = chains, family = Beta(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/prop.litter.rem_model", file_refit = "on_change")
# dat$days.fc <- as.numeric(as.factor(dat$days))
#
# # monotonic effect of time
# fit_mo <- brm(prop.litter.rem ~ trat * mo(days.fc) + (1 | plot.f), data = dat, chains = chains, family = Beta(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/prop.litter.rem_model_monotonic", file_refit = "on_change")
```
```{r, results='asis'}
extended_summary(read.file = "./data/processed/prop.litter.rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
## Remaining wood
```{r, eval = TRUE, out.width = "100%", echo = FALSE, fig.align= "center"}
agg_rem_w <- aggregate(perc.wood.rem ~ colecta + trat + days, dat, mean)
agg_rem_w$sd <- aggregate(perc.wood.rem ~ colecta + trat + days, dat, sd)$perc.wood.rem
agg_rem_w$se <- aggregate(perc.wood.rem ~ colecta + trat + days, dat, se)$perc.wood.rem
pd <- position_dodge(15)
ggplot(agg_rem_w, aes(x = days, y = perc.wood.rem, color = trat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc.wood.rem + se, ymin = perc.wood.rem - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Wood mass remaining (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_rem_w$days),
labels = unique(agg_rem_w$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/wood_mass_remaining.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/wood_mass_remaining.tiff)
```{r, eval = FALSE}
dat$prop.wood.rem <- dat$perc.wood.rem / 100
fit2 <- brm(prop.wood.rem ~ trat * days.sc + (1 | plot.f), data = dat, chains = chains, family = Beta(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/prop.wood.rem_model", file_refit = "on_change")
```
```{r, results='asis'}
extended_summary(read.file = "./data/processed/prop.wood.rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
# K
## Litter
Correlation between Silvia's and Andrea's K
```{r, eval = TRUE}
k_vals <- read.csv("./data/raw/k-values-corr.csv")
cor(k_vals$av.k.litter, k_vals$sil.k.litter)
cor(k_vals$av.k.wood, k_vals$sil.k.wood)
```
Comparing all treatments vs control
```{r}
agg_kvals <- aggregate(sil.k.litter ~ Treatment, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.litter ~ Treatment, k_vals, se)$sil.k.litter
agg_kvals$sd <- aggregate(sil.k.litter ~ Treatment, k_vals, sd)$sil.k.litter
agg_kvals$n <- aggregate(sil.k.litter ~ Treatment, k_vals, length)$sil.k.litter
agg_kvals$treat <- factor(agg_kvals$Treatment, levels = c("C", "N", "P", "NP"))
agg_kvals$n.labels <- paste("n =", agg_kvals$n)
# composed box plot
ggplot(k_vals, aes(x = Treatment, y = sil.k.litter)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the right
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
gghalves::geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)") +
# ylim(c(-0.39, 0.145)) +
geom_text(data = agg_kvals, aes(y = rep(0.4, nrow(agg_kvals)), x = Treatment, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("C" = "Control", "N" = "+N",
"NP" = "+NP", "P" = "+P"))
ggsave("./output/litter_decomposition_k_by_treatment.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/litter_decomposition_k_by_treatment.tiff)
```{r, eval = FALSE}
fit_k.litter <- brm(sil.k.litter ~ Treatment + (1 | quadrat), data = k_vals, chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_litter_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE, read.file = "./data/processed/k_litter_model.rds")
```
Phosphorus vs no-phosphorus
```{r, eval = FALSE}
fit_k.litter.p.np <- brm(sil.k.litter ~ p.treat + (1 | quadrat), data = k_vals, chains = chains, family = gaussian(), iter = iters, , control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_litter_p_nop_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/k_litter_p_nop_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
```{r}
agg_kvals <- aggregate(sil.k.litter ~ p.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.litter ~ p.treat, k_vals, se)$sil.k.litter
agg_kvals$sd <- aggregate(sil.k.litter ~ p.treat, k_vals, sd)$sil.k.litter
agg_kvals$p.treat <- factor(agg_kvals$p.treat, labels = c("No P", "P"))
agg_kvals$n <- aggregate(sil.k.litter ~ p.treat, k_vals, length)$sil.k.litter
agg_kvals$n.labels <- paste("n =", agg_kvals$n)
k_vals$p.treat <- factor(k_vals$p.treat, labels = c("No P", "P"))
# composed box plot
rn_p_litter <- ggplot(k_vals, aes(x = p.treat, y = sil.k.litter)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
gghalves::geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)", x = "P treatment") +
geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = p.treat, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))
rn_p_litter
ggsave("./output/litter_phosphorus_decomposition_k.tiff", dpi = 300)
```
Nitrogen vs no-nitrogen
```{r, eval = FALSE}
fit_k.litter.n.nn <- brm(sil.k.litter ~ n.treat + (1 | quadrat), data = k_vals, chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_litter_n_no_n_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/k_litter_n_no_n_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
```{r}
agg_kvals <- aggregate(sil.k.litter ~ n.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.litter ~ n.treat, k_vals, se)$sil.k.litter
agg_kvals$sd <- aggregate(sil.k.litter ~ n.treat, k_vals, sd)$sil.k.litter
agg_kvals$n.treat <- factor(agg_kvals$n.treat, labels = c("No N", "N"))
agg_kvals$n <- aggregate(sil.k.litter ~ n.treat, k_vals, length)$sil.k.litter
agg_kvals$n.labels <- paste("n =", agg_kvals$n)
k_vals$n.treat <- factor(k_vals$n.treat, labels = c("No N", "N"))
# composed box plot
rn_n_litter <- ggplot(k_vals, aes(x = n.treat, y = sil.k.litter)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
gghalves::geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)", x = "P treatment") +
geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))
rn_n_litter
ggsave("./output/litter_nitrogen_decomposition_k.tiff", dpi = 300)
```
## Wood
Comparing all treatments vs control
```{r}
agg_kvals <- aggregate(sil.k.wood ~ Treatment, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.wood ~ Treatment, k_vals, se)$sil.k.wood
agg_kvals$sd <- aggregate(sil.k.wood ~ Treatment, k_vals, sd)$sil.k.wood
agg_kvals$n <- aggregate(sil.k.wood ~ Treatment, k_vals, length)$sil.k.wood
agg_kvals$treat <- factor(agg_kvals$Treatment, levels = c("C", "N", "P", "NP"))
agg_kvals$n.labels <- paste("n =", agg_kvals$n)
# composed box plot
ggplot(k_vals, aes(x = Treatment, y = sil.k.wood)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
gghalves::geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)") +
# ylim(c(-0.39, 0.145)) +
geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = Treatment, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("C" = "Control", "N" = "+N",
"NP" = "+NP", "P" = "+P"))
ggsave("./output/wood_decomposition_k_by_treatment.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/wood_decomposition_k_by_treatment.tiff)
```{r, eval = FALSE}
fit_k.wood.p.treat <- brm(sil.k.wood ~ Treatment + (1 | quadrat), data = k_vals, chains = chains, family = gaussian(), iter = iters * 2, control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_wood_treatment_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/k_wood_treatment_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
# Phosphorus vs no-phosphorus
```{r}
agg_kvals <- aggregate(sil.k.wood ~ p.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.wood ~ p.treat, k_vals, se)$sil.k.wood
agg_kvals$sd <- aggregate(sil.k.wood ~ p.treat, k_vals, sd)$sil.k.wood
agg_kvals$p.treat <- factor(agg_kvals$p.treat, labels = c("No P", "P"))
agg_kvals$n <- aggregate(sil.k.wood ~ p.treat, k_vals, length)$sil.k.wood
agg_kvals$n.labels <- paste("n =", agg_kvals$n)
k_vals$p.treat <- factor(k_vals$p.treat, labels = c("No P", "P"))
# composed box plot
rn_p_wood <- ggplot(k_vals, aes(x = p.treat, y = sil.k.wood)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
gghalves::geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)", x = "P treatment") +
geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = p.treat, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))
rn_p_wood
ggsave("./output/phosphorus_decomposition_k.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/phosphorus_decomposition_k.tiff)
```{r, eval = FALSE}
fit_k.wood.p.no.p <- brm(sil.k.wood ~ p.treat + (1 | quadrat), data = k_vals, chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_wood_p_no_p_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/k_wood_p_no_p_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
# Nitrogen vs no-nitrogen
```{r}
agg_kvals <- aggregate(sil.k.wood ~ n.treat, k_vals, mean)
agg_kvals$se <- aggregate(sil.k.wood ~ n.treat, k_vals, se)$sil.k.wood
agg_kvals$sd <- aggregate(sil.k.wood ~ n.treat, k_vals, sd)$sil.k.wood
agg_kvals$n.treat <- factor(agg_kvals$n.treat, labels = c("No N", "N"))
agg_kvals$n <- aggregate(sil.k.wood ~ n.treat, k_vals, length)$sil.k.wood
agg_kvals$n.labels <- paste("n =", agg_kvals$n)
k_vals$n.treat <- factor(k_vals$n.treat, labels = c("No N", "N"))
# composed box plot
rn_n_wood <- ggplot(k_vals, aes(x = n.treat, y = sil.k.wood)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
gghalves::geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)", x = "P treatment") +
geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("No P" = "No P", "P" = "+P"))
rn_n_wood
ggsave("./output/nitrogen_decomposition_k.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/nitrogen_decomposition_k.tiff)
```{r, eval = FALSE}
fit_k.wood.n.no.n <- brm(sil.k.wood ~ n.treat + (1 | quadrat), data = k_vals, chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/k_wood_n_no_n_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/k_wood_n_no_n_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
# Nutrient content by remaining litter mass
## Nitrogen
```{r}
nutr$prom.hoja.reman.sc <- scale(nutr$prom.hoja.reman)
ggplot(data = nutr, aes(x = prom.hoja.reman, y = perc.n, color = treat)) +
geom_point() + labs(x = "Remaining litter mass (% initial)", y = "Litter N concentration (%)", color = "Treatment") +
scale_color_viridis_d(alpha = 0.5) +
geom_smooth(method = "lm", se = FALSE) +
scale_x_reverse()
ggsave("./output/nitrogen_by_litter_mass.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/nitrogen_by_litter_mass.tiff)
```{r, eval = FALSE}
fit_pern_by_rem_leave <- brm(perc.n ~ prom.hoja.reman.sc * treat + (1 | plot), data = nutr, chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/n_per_by_leave_rem_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/n_per_by_leave_rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
## Phosphorus
```{r}
ggplot(data = sub.nutr, aes(x = prom.hoja.reman, y = perc.p, color = treat)) +
geom_point() + labs(x = "Remaining litter mass (% initial)", y = "Litter P concentration (%)", color = "Treatment") +
scale_color_viridis_d(alpha = 0.5) +
geom_smooth(method = "lm", se = FALSE) +
scale_x_reverse()
ggsave("./output/phosphorus_by_litter_mass.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/phosphorus_by_litter_mass.tiff)
```{r, eval = FALSE}
sub.nutr$prom.hoja.reman.sc <- scale(sub.nutr$prom.hoja.reman)
fit_perp_by_rem_leave <- brm(perc.p ~ prom.hoja.reman.sc * treat + (1 | plot), data = sub.nutr, chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/p_per_by_leave_rem_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/p_per_by_leave_rem_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
<!-- light brown box -->
<div class="alert alert-success">
# Takeaways {.unnumbered .unlisted}
- Nitrogen percentage, but no Phosphorus percentage, increases along with remaining litter mass
</div>
# Nutrient ratios
## C:N
```{r}
# excluding plot 9
agg_cn <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, mean)
agg_cn$sd <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, sd)$cn.mol.kg
agg_cn$se <- aggregate(cn.mol.kg ~ colecta + treat + days, nutr, se)$cn.mol.kg
agg_cn$treat <- factor(agg_cn$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_cn, aes(x = days, y = cn.mol.kg, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = cn.mol.kg + se, ymin = cn.mol.kg - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter C:N ratio", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_cn$days),
labels = unique(agg_cn$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/cn_ratio_through_time.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/cn_ratio_through_time.tiff)
```{r, eval = FALSE}
fit.cn <- brm(cn.mol.kg ~ treat * days.sc + (1 | plot.f), data = nutr, chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/cn_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/cn_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
## N:P
```{r}
# excluding plot 9
agg_np <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr, mean)
agg_np$sd <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr, sd)$np.mol.kg
agg_np$se <- aggregate(np.mol.kg ~ colecta + treat + days, sub.nutr, se)$np.mol.kg
agg_np$treat <- factor(agg_np$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_np, aes(x = days, y = np.mol.kg, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = np.mol.kg + se, ymin = np.mol.kg - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter N:P ratio", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_np$days),
labels = unique(agg_np$days)) + theme(legend.position = c(0.2, 0.8))
ggsave("./output/np_ratio_through_time.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/np_ratio_through_time.tiff)
```{r, eval = FALSE}
fit.np <- brm(np.mol.kg ~ treat * days.sc + (1 | plot.f), data = sub.nutr, chains = chains, family = gaussian(), iter = iters, control = list(adapt_delta=0.99, max_treedepth=15), cores = chains, prior = prior, file = "./data/processed/np_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/np_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
## C:P
```{r}
# excluding plot 9
agg_cp <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr, mean)
agg_cp$sd <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr, sd)$cp.mol.kg
agg_cp$se <- aggregate(cp.mol.kg ~ colecta + treat + days, sub.nutr, se)$cp.mol.kg
agg_cp$treat <- factor(agg_cp$treat, levels = c("C", "N", "P", "NP"))
pd <- position_dodge(15)
ggplot(agg_cp, aes(x = days, y = cp.mol.kg, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = cp.mol.kg + se, ymin = cp.mol.kg - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5) +
labs(x= "Time (days)", y = "Litter C:P ratio", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_cp$days),
labels = unique(agg_cp$days)) + theme(legend.position = c(0.9, 0.8))
ggsave("./output/cp_ratio_through_time.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/cp_ratio_through_time.tiff)
```{r, eval = FALSE}
fit.cp <- brm(cp.mol.kg ~ treat * days.sc + (1 | plot.f), data = sub.nutr, chains = chains, family = gaussian(), iter = iters, cores = chains, prior = prior, file = "./data/processed/cp_model", file_refit = "on_change")
```
```{r, eval = TRUE, results='asis'}
extended_summary(read.file = "./data/processed/cp_model.rds", gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE, remove.intercepts = TRUE)
```
---
# Combined plots
## Remaining mass
```{r, out.width="100%", out.height="140%"}
agg_rem_w$substrate <- "Wood"
agg_rem$substrate <- "Litter"
names(agg_rem) <- names(agg_rem_w) <- c("colecta", "trat", "days", "perc.rem", "sd", "se", "substrate")
agg_rem_pooled <- rbind(agg_rem, agg_rem_w)
ggplot(agg_rem_pooled, aes(x = days, y = perc.rem, color = trat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc.rem + se, ymin = perc.rem - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N", "+P", "+NP")) +
labs(x= "Time (days)", y = "Remaining mass (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_rem_pooled$days),
labels = unique(agg_rem_pooled$days)) +
facet_wrap(~ substrate, nrow = 2)
ggsave("./output/remaining_mass_combined.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/remaining_mass_combined.tiff)
## Litter content
```{r, out.width="100%", out.height="140%"}
agg_p$nutrient <- "P"
agg_n$nutrient <- "N"
names(agg_n) <- names(agg_p) <- c("colecta", "treat", "days", "perc", "sd", "se", "nutrient")
agg_np2 <- rbind(agg_n, agg_p)
ggplot(agg_np2, aes(x = days, y = perc, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc + se, ymin = perc - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N", "+P", "+NP")) +
labs(x= "Time (days)", y = "Litter nutrient content (% initial)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_p$days),
labels = unique(agg_p$days)) +
facet_wrap(~ nutrient, nrow = 2)
ggsave("./output/litter_content_combined.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/litter_content_combined.tiff)
## Concentration
```{r, out.width="100%", out.height="140%"}
agg_conc_p$nutrient <- "P"
agg_conc_n$nutrient <- "N"
names(agg_conc_n) <- names(agg_conc_p) <- c("colecta", "treat", "days", "perc", "sd", "se", "nutrient")
agg_conc_np <- rbind(agg_conc_n, agg_conc_p)
ggplot(agg_conc_np, aes(x = days, y = perc, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = perc + se, ymin = perc - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N", "+P", "+NP")) +
labs(x= "Time (days)", y = "Litter nutrient concentration (%)", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_conc_p$days),
labels = unique(agg_conc_p$days)) +
facet_wrap(~ nutrient, nrow = 2, scales = "free_y")
ggsave("./output/concentration_combined.tiff", dpi = 300)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/concentration_combined.tiff)
## Nutrient ratios
```{r, out.width="100%", fig.height=12}
agg_cn$nutrients <- "C:N"
agg_cp$nutrients <- "C:P"
agg_np$nutrients <- "N:P"
names(agg_cn) <- names(agg_cp) <- names(agg_np) <- c("colecta", "treat", "days", "mol.kg", "sd", "se", "nutrients")
agg_ratios <- rbind(agg_cn[, names(agg_cn)], agg_cp[, names(agg_cn)], agg_np[, names(agg_cn)])
ggplot(agg_ratios, aes(x = days, y = mol.kg, color = treat)) +
geom_point(size= 2, position = pd) +
geom_errorbar(aes(ymax = mol.kg + se, ymin = mol.kg - se), width = 0, position = pd) +
geom_line(size = 1.2, position = pd) +
scale_color_viridis_d(alpha = 0.5, labels = c("Control", "+N", "+P", "+NP")) +
labs(x= "Time (days)", y = "Litter nutrient ratios", color = "Treatment") +
scale_x_continuous(breaks = unique(agg_ratios$days),
labels = unique(agg_ratios$days)) +
facet_wrap(~ nutrients, nrow = 3, scales = "free_y")
ggsave("./output/nutrient_ratios_combined.tiff", dpi = 300, width = 10, height = 12)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/nutrient_ratios_combined.tiff)
## Decomposition constant
```{r, out.width="100%", fig.height=12}
rn_n_wood_dat <- rn_n_wood$data[, c("n.treat", "sil.k.wood")]
rn_n_litter_dat <- rn_n_litter$data[, c("n.treat", "sil.k.litter")]
rn_p_wood_dat <- rn_p_wood$data[, c("p.treat", "sil.k.wood")]
rn_p_litter_dat <- rn_p_litter$data[, c("p.treat", "sil.k.litter")]
names(rn_p_wood_dat) <- names(rn_n_wood_dat) <- names(rn_p_litter_dat) <- names(rn_n_litter_dat) <- c("treatment", "k")
rn_p_wood_dat$substrate <- rn_n_wood_dat$substrate <- "wood"
rn_p_litter_dat$substrate <- rn_n_litter_dat$substrate <- "litter"
rn_n_wood_dat$nutrient <- rn_n_litter_dat$nutrient <- "N"
rn_p_litter_dat$nutrient <- rn_p_wood_dat$nutrient <- "P"
rn_k_dat <- rbind(rn_n_wood_dat, rn_p_wood_dat, rn_n_litter_dat, rn_p_litter_dat)
# composed box plot
ggplot(rn_k_dat, aes(x = treatment, y = k)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
gghalves::geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)", x = "P treatment") +
# geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("No P" = "No P", "P" = "+P")) +
facet_grid(~ substrate)
ggsave("./output/decomposition_k_combined_v1.tiff", dpi = 300, width = 10, height = 12)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/decomposition_k_combined_v1.tiff)
```{r}
# composed box plot
ggplot(rn_k_dat, aes(x = treatment, y = k)) +
## add half-violin from {ggdist} package
ggdist::stat_halfeye(
fill = fill_color,
alpha = 0.5,
## custom bandwidth
adjust = .5,
## adjust height
width = .6,
.width = 0,
## move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(fill = fill_color,
width = .15,
## remove outliers
outlier.shape = NA ## `outlier.shape = NA` works as well
) +
## add justified jitter from the {gghalves} package
gghalves::geom_half_point(
color = fill_color,
## draw jitter on the left
side = "l",
## control range of jitter
range_scale = .4,
## add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
labs(y = "Decomposition constant (k)", x = "P treatment") +
# geom_text(data = agg_kvals, aes(y = rep(-0.387, nrow(agg_kvals)), x = n.treat, label = n.labels), nudge_x = 0, size = 6) +
theme_classic(base_size = 18) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels=c("No P" = "No P", "P" = "+P")) +
facet_grid(substrate ~ nutrient, scales = "free_x")
ggsave("./output/decomposition_k_combined_v12.tiff", dpi = 300, width = 10, height = 12)
```
[Download image](https://github.com/maRce10/litter_decomposition_EFFEX/raw/master/output/decomposition_k_combined_v2.tiff)
# Combined model diagnostics
```{r, eval = TRUE}
check_rds_fits(path = "./data/processed", html = TRUE)
```
<!-- add packages used, system details and versions -->
# Session information {.unnumbered .unlisted}
```{r session info, echo=F}
sessionInfo()
```